clear
set type double
version 16.1


/* firm */
*************** false discovery rate *** added on 2may23 * (include firms with g28_plus ==.) ****************
global alpha = .03
use ../firmleveltest/testresults_020119.dta, clear
keep if regionfirmTG == 1
*drop  id rank1* dep* indep* tie* rank_* fy_con
compress
count if (( ftt_ts1s > 1.64 & ftt_ts1s < .)) 
count if ((ftt_ts1p < -1.64 & ftt_ts1p < .)) 
gen firm1143 = ((ftt_ts1p < -1.64 & ftt_ts1p < . )|( ftt_ts1s > 1.64 & ftt_ts1s < .)) 
count if firm1143 == 1
gen normal01 = invnormal(runiform()) if regionfirmTG == 1


foreach x in "ftt_ts1s" "ftt_ts1p"{
	preserve
		*keep if regionfirmTG == 1 & ftt_N1 < .
		*local x = "ftt_ts1s"
		*local x = "ftt_ts1p"
		keep if regionfirmTG == 1 & ftt_N1 > 10 & ftt_N1 < . 
		count
		gen pv = 1-normal(`x')
		sum partic ftt_N1 pv

		hist pv , bin(100)
		*apply shrinkage 15may23
		gen mod_pv = pv
		/*replace mod_pv = pv + 1 if ftt_tau1s < 0
		replace mod_pv = mod_pv/2 
		cap graph drop gmodpv*/
		*hist mod_pv, bin(100) name(gmodpv)
		gen Yb = `x'
		gen mpv = 1 - pv
		*egen sigma = std(Yb)';
		gen Y1 = Yb
		egen meanY1 = mean(Y1)
		egen A = mean((Y1-meanY1)^2 - 1)
		count if Y1 < .
		gen B = (2*(1/r(N))^2 * r(N))
		egen mu2 = rowmax(A B)

		gen shrinkY = meanY1 + (mu2/(mu2+1)) * (Y1-meanY1)
		gen shrinkpv = 1-normal(shrinkY)
		twoway(hist pv, bin(100) color(navy))(hist shrinkpv, bin(100) color(%40)), legend(off)
		graph export "../graph/fdr_hist_`x'.png", as(png) replace 

		sort `x'
		count if `x' < .
		local N = r(N)
		cap drop n 
		qui gen n = _n / `N' if `x' < .
		cap drop mshrinkpv
		gen mshrinkpv = 1-shrinkpv
		*gen normal = invnormal(_n)
		*gen cdf_shrinkY = / r(N)
		if "`x'" == "ftt_ts1p"{
			count if shrinkY < -1.65
		}
		else{
			count if shrinkY < 1.65
		}
		local reject95 = r(N) / `N'
		di `reject95'	
		
		***** shinkage graph *****
		if "`x'" == "ftt_ts1p"{
			local str_reject95 = string(`reject95', "%6.2f")
			twoway(line mshrinkpv shrinkY, color(gs13) lwidth(vthick) xax(1) text(1.06 -1.7 "-1.65", size(small))) ///
			(line n shrinkY, color(gs2) lwidth(vthick) ///
					text(.05 -4.5 "0.05", size(small)) text(`reject95' -4.5 "`str_reject95'", size(small)) ) ///
				if abs(shrinkY) < 4, graphregion(color(white)) legend(off) xtitle("") ///
				title("") ylabel(,angle(0)) xlabel(-4(2)4) xline(-4 -2 0 2 4, lcolor(gs13) ///
				lpattern(solid) lwidth(vthin) ax(1)) /*xlabel("",ax(2)) xtitle("", ax(2))*/ ///
				xline(-1.65, lcolor(gs2) lpattern(solid) lwidth(medium) ax(1) lstyle(foreground))  ///
				xline(0, lcolor(gs13) lpattern(solid) lwidth(thin) ax(1)) ///
				yline(.05, lcolor(gs8) lpattern(shortdash) lwidth(vthin) ) ///
				yline(`reject95', lcolor(gs8) lpattern(shortdash) lwidth(vthin)) ///
				ylabel(0(.5)1) scale(1.2)
		}
		else{
			twoway(line mshrinkpv shrinkY, color(gs13) lwidth(vthick) xax(1) text(1.06 1.65 "1.65", size(small))) ///
			(line n shrinkY, color(gs2) lwidth(vthick) ///
					text(.94 -4.5 "0.95", size(small)) text(.85 -4.5 "0.86", size(small)) ) ///
				if abs(shrinkY) < 4, graphregion(color(white)) legend(off) xtitle("") ///
				title("") ylabel(,angle(0)) xlabel(-4(2)4) xline(-4 -2 0 2 4, lcolor(gs13) ///
				lpattern(solid) lwidth(vthin) ax(1)) /*xlabel("",ax(2)) xtitle("", ax(2))*/ ///
				xline(1.65, lcolor(gs2) lpattern(solid) lwidth(medium) ax(1) lstyle(foreground))  ///
				xline(0, lcolor(gs13) lpattern(solid) lwidth(thin) ax(1)) ///
				yline(.95, lcolor(gs8) lpattern(shortdash) lwidth(vthin) ) ///
				yline(.859, lcolor(gs8) lpattern(shortdash) lwidth(vthin)) ///
				ylabel(0(.5)1) scale(1.2)
		}
		sleep 10
			
		*no shrinkage (added 8/18/2023)
		if "`x'" == "ftt_ts1p"{
			count if Y1 < -1.65 
			local reject95 = r(N) / `N'
			local str_reject95 = string(`reject95', "%6.2f")
			di `reject95'
				twoway(line mpv Y1, color(gs13) lwidth(vthick) xax(1) text(1.06 -1.7 "-1.65", size(small))) ///
			(line n Y1, color(gs2) lwidth(vthick) ///
					text(.05 -4.5 "0.05", size(small)) text(`reject95' -4.5 "`str_reject95'", size(small)) ) ///
				if abs(Y1) < 4, graphregion(color(white)) legend(off) xtitle("") ///
				title("") ylabel(,angle(0)) xlabel(-4(2)4) xline(-4 -2 0 2 4, lcolor(gs13) ///
				lpattern(solid) lwidth(vthin) ax(1)) /*xlabel("",ax(2)) xtitle("", ax(2))*/ ///
				xline(-1.65, lcolor(gs2) lpattern(solid) lwidth(medium) ax(1) lstyle(foreground))  ///
				xline(0, lcolor(gs13) lpattern(solid) lwidth(thin) ax(1)) ///
				yline(.05, lcolor(gs8) lpattern(shortdash) lwidth(vthin) ) ///
				yline(`reject95', lcolor(gs8) lpattern(shortdash) lwidth(vthin)) ///
				ylabel(0(.5)1) scale(1.2)		
		}
		else{
			count if Y1 < 1.65
			local reject95 = r(N) / `N'
			local str_reject95 = string(`reject95', "%6.2f")
			di `reject95'
				twoway(line mpv Y1, color(gs13) lwidth(vthick) xax(1) text(1.06 1.65 "1.65", size(small))) ///
			(line n Y1, color(gs2) lwidth(vthick) ///
					text(.94 -4.5 "0.95", size(small)) text(`reject95' -4.5 "`str_reject95'", size(small)) ) ///
				if abs(Y1) < 4, graphregion(color(white)) legend(off) xtitle("") ///
				title("") ylabel(,angle(0)) xlabel(-4(2)4) xline(-4 -2 0 2 4, lcolor(gs13) ///
				lpattern(solid) lwidth(vthin) ax(1)) /*xlabel("",ax(2)) xtitle("", ax(2))*/ ///
				xline(1.65, lcolor(gs2) lpattern(solid) lwidth(medium) ax(1) lstyle(foreground))  ///
				xline(0, lcolor(gs13) lpattern(solid) lwidth(thin) ax(1)) ///
				yline(.95, lcolor(gs8) lpattern(shortdash) lwidth(vthin) ) ///
				yline(`reject95', lcolor(gs8) lpattern(shortdash) lwidth(vthin)) ///
				ylabel(0(.5)1) scale(1.2)
			

		}
		graph export "../graph/fdr_noshrinkage_cdfplot_`x'.png", as(png) replace 
		
		
		* computation of pFDR (Storey 2002)
		if "`x'" == "ftt_ts1s"{
			local lambda = -1.65 /*lambda is a tuning parameter */
			count if Y1<-1.65		
			local W = r(N)
			local hat_pi = `W'/(0.05*`N')
			di `hat_pi'
			local hat_pFDR = (`hat_pi'*0.05)/(1-`reject95') 
			di `hat_pFDR'
		}
		if "`x'" == "ftt_ts1p"{
			replace mod_pv = 1-mod_pv
		}
		sum mod_pv
		local N = r(N)
		di "N: `N'"
		sort mod_pv
		gen bhstep = (_n / `N') * $alpha
		order regionfirmTG regionfirmID pv mod_pv bhstep
		gsort -mod_pv
		gen reject = .
		replace reject = 1 if ((mod_pv < bhstep) | (reject[_n-1] == 1))
		replace reject = . if mod_pv == .
		order regionfirmTG regionfirmID pv mod_pv bhstep reject
		sort mod_pv
		gen bkystep = (_n / `N') * ($alpha / (1 + $alpha)) 
		gen rejectbky = .
		gsort -mod_pv	
		replace rejectbky = 1 if ((mod_pv < bkystep) | (reject[_n-1] == 1))
		replace rejectbky = . if mod_pv == .
		order regionfirmTG regionfirmID mod_pv pv reject bkystep rejectbky
		count if rejectbky== 1
		local N_bky = r(N)
		sort mod_pv	
		gen bkystep_final = (_n / (`N' - `N_bky')) * ($alpha) 
		gsort -mod_pv
		gen rejectbky_final = .
		replace rejectbky_final = 1 if ((mod_pv < bkystep_final) | (reject[_n-1] == 1))
		replace rejectbky_final = . if mod_pv == .
		count if reject == 1
		count if rejectbky_final == 1
	restore
}
